

#---Portmann & Stojanovic 2017: Electoral discrimination against immigrant-origin candidates---#

#load("PrefVotes")

#-------------------- ESTIMATE MODELS (BETA DISTRIBUTION, LOGIT LINK)

#Models 

model1 <- glmmadmb(crossoff_pct ~ nonswissnames + (~1|list_id), data=d, family="Beta")

model2 <- glmmadmb(crossoff_pct ~ nonswissnames + incumbent + c_relposition_rc + precumulated + 
                     sex_d + age_c3 + CodeBerufK3 + num_cand + party_c3 + commune
                     (~1|list_id), data=d, family="Beta", admb.opts=admbControl(shess=FALSE,noinit=FALSE))

model3 <- glmmadmb(crossoff_pct ~ nonswissnames + incumbent + c_relposition_rc + precumulated + 
                     sex_d + age_c3 + CodeBerufK3 + num_cand + party_c3 + nonswissnames:party_c3 + commune + 
                   (~1|list_id), data=d, family="Beta", admb.opts=admbControl(shess=FALSE,noinit=FALSE))

#-------------------- COEFFICIENT PLOT

#Put data from models in temp. data frame

model1Frame <- data.frame(Variable = rownames(summary(model1)$coef),
                          Coefficient = summary(model1)$coef[, 1],
                          SE = summary(model1)$coef[, 2],
                          Model = "Model 1")
model2Frame <- data.frame(Variable = rownames(summary(model2)$coef),
                          Coefficient = summary(model2)$coef[, 1],
                          SE = summary(model2)$coef[, 2],
                          Model = "Model 2")
model3Frame <- data.frame(Variable = rownames(summary(model3)$coef),
                          Coefficient = summary(model3)$coef[, 1],
                          SE = summary(model3)$coef[, 2],
                          Model = "Model 3")

#Combine data.frames

allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))

#Specify 95% confidence intervals

interval1 <- -qnorm((1-0.95)/2)  # 95% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

#Change order of variables

order_var <- c("(Intercept)", "nonswissnames1:party_c32", 
               "nonswissnames1:party_c31", "party_c32", "party_c31", "num_cand", "CodeBerufK33", "CodeBerufK32", "age_c32", "age_c31", "sex_d1", 
               "precumulated1","c_relposition_rc","incumbent1","nonswissnames1")

allModelFrame$Variable <- factor(allModelFrame$Variable, levels=order_var)

#Drop commune fixed effects

cc=is.na(allModelFrame$Variable)
m=which(cc==c("TRUE"))
allModelFrame=allModelFrame [-m,]

#Define labels

my_labels = c("(Intercept)" = "Intercept", "nonswissnames1:party_c32" = "Non-Swiss name*Party position = right", 
              "nonswissnames1:party_c31" = "Non-Swiss name*Party position = centre-right",
              "num_cand" = "Number of candidates on ballot","party_c32"="Party position: right", "party_c31"="Party position: centre-right", 
              "nonswissnames1"="Non-Swiss name", "age_c32"="Age: 50+ years", "age_c31"="Age: 30-50 years", "sex_d1"= "Female", 
              "precumulated1" = "Pre-cumulated","c_relposition_rc" = "Relative position on ballot",
              "CodeBerufK32" = "Profession: medium skilled", "CodeBerufK33" = "Profession: low skilled",
              "incumbent1" = "Incumbent")

#Define theme

my_theme <- theme(axis.text.y = element_text(size = 12, colour = "black")) + theme(axis.text.x= element_text(size = 12, colour = "black")) +  
  theme(axis.title = element_text(size = 12)) + theme(legend.position="bottom") + theme(legend.title = element_text(size = 12)) + 
  theme(legend.text = element_text(size = 12)) + theme(plot.margin = unit(c(0,1,0,0), "cm")) 

#Coefficient plot

z <- PlotCoefficients(z)
z <- PlotCoefficients(z) + scale_y_continuous("Effect on negative preference votes", limits = c(-2, 2)) + scale_x_discrete("", labels = my_labels) +
  my_theme + theme(axis.text.y = element_text(face = c("plain", "bold", "bold", "plain", "plain", "plain", "plain", "plain", "plain", "plain", "plain", "plain", "plain", "plain", "bold"))) +
  scale_colour_manual("Model", values = c("gray60", "gray40", "black")) + theme(legend.title = element_blank())
z

tiff("Figure 1.tiff", width = 8.4, height = 7, units = 'in', res = 300)
plot(z)
dev.off()



